"""
DaoMath Astrophysics (道数天体物理)
从 ±0 先天均衡构造宇宙理论

天体物理的 ±0 诠释:
    引力吸引 ↔ 宇宙膨胀 (对偶)
    物质 ↔ 暗物质 (对偶)
    能量 ↔ 暗能量 (对偶)
    恒星诞生 ↔ 恒星死亡 (对偶)
    
理论层级:
    ±0 (先天均衡)
     ↓
    时空几何 (广义相对论)
     ↓
    恒星结构与演化
     ↓
    星系动力学
     ↓
    宇宙学 (膨胀, 暗物质, 暗能量)
     ↓
    大爆炸理论
     ↓
    结构形成
"""

from __future__ import annotations
from typing import List, Tuple, Callable, Optional, Dict, Any, Union
from dataclasses import dataclass, field
from abc import ABC, abstractmethod
from enum import Enum
import numpy as np
import scipy.linalg as la
from scipy.integrate import odeint, solve_ivp, quad
from scipy.optimize import fsolve, brentq
from scipy.special import spherical_jn, legendre
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from core.axioms import InnateBalance, DualPair


# ============================================================================
# 物理常数（天文单位）
# ============================================================================

class AstroConstants:
    """天体物理常数"""
    # 基本常数
    C = 2.99792458e8                    # 光速 (m/s)
    G = 6.67430e-11                     # 引力常数 (m³/kg/s²)
    H_BAR = 1.054571817e-34             # ℏ (J·s)
    K_B = 1.380649e-23                  # Boltzmann 常数 (J/K)
    SIGMA_SB = 5.670374419e-8           # Stefan-Boltzmann 常数 (W/m²/K⁴)
    
    # 天文单位
    AU = 1.495978707e11                 # 天文单位 (m)
    PARSEC = 3.0857e16                  # 秒差距 (m)
    LIGHTYEAR = 9.4607e15               # 光年 (m)
    
    # 太阳参数
    M_SUN = 1.98892e30                  # 太阳质量 (kg)
    R_SUN = 6.96e8                      # 太阳半径 (m)
    L_SUN = 3.828e26                    # 太阳光度 (W)
    T_SUN = 5778                        # 太阳表面温度 (K)
    
    # 地球参数
    M_EARTH = 5.972e24                  # 地球质量 (kg)
    R_EARTH = 6.371e6                   # 地球半径 (m)
    
    # 宇宙学常数
    H0 = 67.4                           # Hubble 常数 (km/s/Mpc)
    OMEGA_M = 0.315                     # 物质密度参数
    OMEGA_LAMBDA = 0.685                # 暗能量密度参数
    OMEGA_B = 0.049                     # 重子物质密度参数
    OMEGA_R = 9.0e-5                    # 辐射密度参数
    T_CMB = 2.725                       # 宇宙微波背景温度 (K)
    
    # 其他
    M_PLANCK = 2.176434e-8              # Planck 质量 (kg)
    L_PLANCK = 1.616255e-35             # Planck 长度 (m)
    T_PLANCK = 5.391247e-44             # Planck 时间 (s)


# ============================================================================
# Level 1: 恒星结构
# ============================================================================

class StellarStructure:
    """
    恒星结构方程
    
    描述恒星内部的平衡状态
    
    基于 ±0 的理解:
    - 引力压缩 ↔ 气体压力（流体静力学平衡）
    - 核反应产能 ↔ 表面辐射（能量平衡）
    """
    
    def __init__(self, 
                 mass: float,
                 radius: float,
                 luminosity: float,
                 composition: Dict[str, float]):
        """
        Args:
            mass: 恒星质量 (M_☉)
            radius: 恒星半径 (R_☉)
            luminosity: 光度 (L_☉)
            composition: 化学组成 {'H': X, 'He': Y, 'Z': Z}
        """
        self.M = mass * AstroConstants.M_SUN
        self.R = radius * AstroConstants.R_SUN
        self.L = luminosity * AstroConstants.L_SUN
        self.X = composition.get('H', 0.70)    # 氢丰度
        self.Y = composition.get('He', 0.28)   # 氦丰度
        self.Z = composition.get('Z', 0.02)    # 金属丰度
    
    def hydrostatic_equilibrium(self, r: float, M_r: float, rho: float) -> float:
        """
        流体静力学平衡
        
        dP/dr = -G M(r) ρ(r) / r²
        
        基于 ±0: 压力梯度(+) ↔ 引力(-)
        
        Args:
            r: 半径 (m)
            M_r: 包含质量 (kg)
            rho: 密度 (kg/m³)
        
        Returns:
            压力梯度 dP/dr (Pa/m)
        """
        G = AstroConstants.G
        
        if r == 0:
            return 0
        
        return -G * M_r * rho / r**2
    
    def mass_continuity(self, r: float, rho: float) -> float:
        """
        质量连续性方程
        
        dM/dr = 4π r² ρ(r)
        
        Args:
            r: 半径 (m)
            rho: 密度 (kg/m³)
        
        Returns:
            质量梯度 dM/dr (kg/m)
        """
        return 4 * np.pi * r**2 * rho
    
    def energy_generation(self, rho: float, T: float) -> float:
        """
        能量产生率（核聚变）
        
        ε ~ ρ T^ν
        
        - PP 链: ν ≈ 4
        - CNO 循环: ν ≈ 18
        
        Args:
            rho: 密度 (kg/m³)
            T: 温度 (K)
        
        Returns:
            能量产生率 ε (W/kg)
        """
        # 主序星（PP 链主导）
        epsilon_0 = 1.0e-5  # W/kg (归一化)
        T_0 = 1.5e7         # K
        
        # 温度依赖（简化）
        if T < 1e7:
            return 0
        
        epsilon = epsilon_0 * rho * (T / T_0)**4
        
        return epsilon
    
    def radiative_energy_transport(self, 
                                   r: float, 
                                   T: float, 
                                   L_r: float,
                                   rho: float,
                                   kappa: float) -> float:
        """
        辐射能量输运
        
        dT/dr = -3 κ ρ L(r) / (16π a c T³ r²)
        
        其中 κ 是不透明度
        
        Args:
            r: 半径 (m)
            T: 温度 (K)
            L_r: 光度 (W)
            rho: 密度 (kg/m³)
            kappa: 不透明度 (m²/kg)
        
        Returns:
            温度梯度 dT/dr (K/m)
        """
        if r == 0:
            return 0
        
        c = AstroConstants.C
        a = 4 * AstroConstants.SIGMA_SB / c  # 辐射常数
        
        dT_dr = -3 * kappa * rho * L_r / (16 * np.pi * a * c * T**3 * r**2)
        
        return dT_dr
    
    def opacity(self, rho: float, T: float) -> float:
        """
        不透明度 κ (Kramers' law)
        
        κ ~ ρ T^(-7/2)
        
        Args:
            rho: 密度 (kg/m³)
            T: 温度 (K)
        
        Returns:
            不透明度 (m²/kg)
        """
        kappa_0 = 0.02  # m²/kg (归一化)
        rho_0 = 1e3     # kg/m³
        T_0 = 1e6       # K
        
        kappa = kappa_0 * (rho / rho_0) * (T / T_0)**(-3.5)
        
        return max(kappa, 1e-4)  # 下限
    
    def solve_structure(self, num_points: int = 100) -> Dict[str, np.ndarray]:
        """
        求解恒星结构（简化数值解）
        
        Returns:
            {'r': 半径, 'M': 质量, 'P': 压力, 'T': 温度, 'rho': 密度}
        """
        r_vals = np.linspace(0.01 * self.R, self.R, num_points)
        
        # 初始猜测（多方球）
        rho_c = 1.5e5  # kg/m³ (中心密度)
        T_c = 1.5e7    # K (中心温度)
        P_c = 2e16     # Pa (中心压力)
        
        M_vals = np.zeros(num_points)
        rho_vals = np.zeros(num_points)
        P_vals = np.zeros(num_points)
        T_vals = np.zeros(num_points)
        
        # 边界条件
        M_vals[0] = 0
        rho_vals[0] = rho_c
        P_vals[0] = P_c
        T_vals[0] = T_c
        
        # 向外积分（简化 Euler 法）
        for i in range(1, num_points):
            r = r_vals[i]
            dr = r_vals[i] - r_vals[i-1]
            
            # 多方关系（简化）
            n = 3  # 多方指数
            K = P_vals[0] / rho_vals[0]**(1 + 1/n)
            
            # 质量
            dM_dr = self.mass_continuity(r, rho_vals[i-1])
            M_vals[i] = M_vals[i-1] + dM_dr * dr
            
            # 压力
            dP_dr = self.hydrostatic_equilibrium(r, M_vals[i], rho_vals[i-1])
            P_vals[i] = P_vals[i-1] + dP_dr * dr
            
            # 密度（多方关系）
            if P_vals[i] > 0:
                rho_vals[i] = (P_vals[i] / K)**(n / (n + 1))
            else:
                rho_vals[i] = 0
            
            # 温度（理想气体）
            mu = 0.6  # 平均分子量
            m_H = 1.67e-27  # kg
            k_B = AstroConstants.K_B
            
            if rho_vals[i] > 0:
                T_vals[i] = P_vals[i] * mu * m_H / (rho_vals[i] * k_B)
            else:
                T_vals[i] = 0
        
        return {
            'r': r_vals,
            'M': M_vals,
            'P': P_vals,
            'T': T_vals,
            'rho': rho_vals
        }
    
    def central_temperature(self) -> float:
        """
        估算中心温度
        
        T_c ~ G M μ m_H / (k_B R)
        
        维里定理
        """
        G = AstroConstants.G
        k_B = AstroConstants.K_B
        mu = 0.6
        m_H = 1.67e-27
        
        T_c = G * self.M * mu * m_H / (k_B * self.R)
        
        return T_c
    
    def surface_temperature(self) -> float:
        """
        表面温度（Stefan-Boltzmann 定律）
        
        L = 4π R² σ T_eff⁴
        """
        sigma = AstroConstants.SIGMA_SB
        
        T_eff = (self.L / (4 * np.pi * self.R**2 * sigma))**0.25
        
        return T_eff


class StellarEvolution:
    """
    恒星演化
    
    主序星 → 红巨星 → （白矮星/中子星/黑洞）
    
    基于 ±0 的理解:
    - 核燃烧(+) ↔ 引力收缩(-)
    - 膨胀 ↔ 收缩
    - 诞生 ↔ 死亡
    """
    
    class EvolutionStage(Enum):
        """演化阶段"""
        PROTOSTAR = "原恒星"
        MAIN_SEQUENCE = "主序星"
        RED_GIANT = "红巨星"
        HORIZONTAL_BRANCH = "水平分支"
        ASYMPTOTIC_GIANT_BRANCH = "渐近巨星分支"
        PLANETARY_NEBULA = "行星状星云"
        WHITE_DWARF = "白矮星"
        SUPERNOVA = "超新星"
        NEUTRON_STAR = "中子星"
        BLACK_HOLE = "黑洞"
    
    @staticmethod
    def main_sequence_lifetime(mass: float) -> float:
        """
        主序寿命
        
        τ_MS ~ M / L ~ M^(-2.5)
        
        质量-光度关系: L ~ M^3.5
        
        Args:
            mass: 质量 (M_☉)
        
        Returns:
            寿命 (年)
        """
        tau_sun = 1e10  # 太阳主序寿命 ~ 100 亿年
        
        tau_MS = tau_sun * mass**(-2.5)
        
        return tau_MS
    
    @staticmethod
    def mass_luminosity_relation(mass: float) -> float:
        """
        质量-光度关系
        
        L ~ M^α
        
        α ≈ 3.5 (主序星)
        
        Args:
            mass: 质量 (M_☉)
        
        Returns:
            光度 (L_☉)
        """
        if mass < 0.43:
            alpha = 2.3
        elif mass < 2:
            alpha = 4.0
        elif mass < 55:
            alpha = 3.5
        else:
            alpha = 1.0
        
        return mass**alpha
    
    @staticmethod
    def chandrasekhar_limit() -> float:
        """
        Chandrasekhar 极限
        
        M_Ch ≈ 1.44 M_☉
        
        白矮星质量上限
        
        基于 ±0: 电子简并压(+) ↔ 引力(-)
        """
        return 1.44  # M_☉
    
    @staticmethod
    def tolman_oppenheimer_volkoff_limit() -> float:
        """
        TOV 极限
        
        M_TOV ≈ 2-3 M_☉
        
        中子星质量上限
        """
        return 2.5  # M_☉ (不确定)
    
    @staticmethod
    def schwarzschild_radius(mass: float) -> float:
        """
        Schwarzschild 半径（黑洞视界）
        
        R_s = 2GM/c²
        
        Args:
            mass: 质量 (M_☉)
        
        Returns:
            Schwarzschild 半径 (m)
        """
        G = AstroConstants.G
        c = AstroConstants.C
        M_sun = AstroConstants.M_SUN
        
        R_s = 2 * G * mass * M_sun / c**2
        
        return R_s
    
    @staticmethod
    def final_state(initial_mass: float) -> EvolutionStage:
        """
        最终状态（根据初始质量）
        
        M < 0.08 M_☉: 褐矮星
        0.08 < M < 0.5 M_☉: 红矮星（长寿）
        0.5 < M < 8 M_☉: 白矮星
        8 < M < 25 M_☉: 中子星
        M > 25 M_☉: 黑洞
        
        Args:
            initial_mass: 初始质量 (M_☉)
        """
        if initial_mass < 8:
            return StellarEvolution.EvolutionStage.WHITE_DWARF
        elif initial_mass < 25:
            return StellarEvolution.EvolutionStage.NEUTRON_STAR
        else:
            return StellarEvolution.EvolutionStage.BLACK_HOLE


class NuclearFusion:
    """
    恒星核聚变
    
    基于 ±0 的理解:
    - 核力吸引(+) ↔ 电磁排斥(-)
    - 质量亏损 → 能量释放
    """
    
    @staticmethod
    def pp_chain() -> Dict[str, Any]:
        """
        质子-质子链
        
        4 ¹H → ⁴He + 2e⁺ + 2νₑ + 2γ
        
        Q = 26.7 MeV
        
        主要发生在 T < 2×10⁷ K
        """
        return {
            'reaction': '4p → He-4 + 2e+ + 2νe + 2γ',
            'Q_value': 26.7,  # MeV
            'temperature_range': (1e7, 2e7),  # K
            'dominant_in': 'M < 1.5 M_☉'
        }
    
    @staticmethod
    def cno_cycle() -> Dict[str, Any]:
        """
        CNO 循环
        
        碳-氮-氧作为催化剂
        
        主要发生在 T > 2×10⁷ K
        """
        return {
            'reaction': '4p + C/N/O → He-4 + C/N/O',
            'Q_value': 25.0,  # MeV
            'temperature_range': (2e7, 1e8),  # K
            'dominant_in': 'M > 1.5 M_☉'
        }
    
    @staticmethod
    def triple_alpha() -> Dict[str, Any]:
        """
        3α 过程
        
        3 ⁴He → ¹²C
        
        发生在红巨星氦燃烧阶段
        """
        return {
            'reaction': '3 He-4 → C-12',
            'Q_value': 7.3,  # MeV
            'temperature_range': (1e8, 2e8),  # K
            'phase': 'Helium burning'
        }
    
    @staticmethod
    def binding_energy_per_nucleon(A: int) -> float:
        """
        每核子结合能
        
        Fe-56 最稳定（峰值）
        
        基于 ±0: 轻核聚变(+能量) ↔ 重核裂变(+能量)
        
        Args:
            A: 质量数
        
        Returns:
            结合能 (MeV/nucleon)
        """
        # 半经验公式（简化）
        if A == 4:
            return 7.07  # He-4
        elif A == 12:
            return 7.68  # C-12
        elif A == 56:
            return 8.79  # Fe-56 (最高)
        else:
            # 近似
            return 8.0 * (1 - abs(A - 56) / 200)


# ============================================================================
# Level 2: 星系动力学
# ============================================================================

class Galaxy:
    """
    星系
    
    恒星、气体、暗物质的引力束缚系统
    
    基于 ±0 的理解:
    - 可见物质 ↔ 暗物质
    - 旋转 ↔ 引力
    """
    
    class GalaxyType(Enum):
        """星系类型（Hubble 分类）"""
        ELLIPTICAL = "椭圆星系"
        SPIRAL = "旋涡星系"
        BARRED_SPIRAL = "棒旋星系"
        IRREGULAR = "不规则星系"
        LENTICULAR = "透镜状星系"
    
    def __init__(self,
                 galaxy_type: GalaxyType,
                 mass: float,
                 radius: float,
                 dark_matter_fraction: float = 0.85):
        """
        Args:
            galaxy_type: 星系类型
            mass: 总质量 (M_☉)
            radius: 特征半径 (kpc)
            dark_matter_fraction: 暗物质比例
        """
        self.type = galaxy_type
        self.M_total = mass * AstroConstants.M_SUN
        self.R = radius * 1e3 * AstroConstants.PARSEC
        self.f_dm = dark_matter_fraction
        
        self.M_visible = self.M_total * (1 - dark_matter_fraction)
        self.M_dark = self.M_total * dark_matter_fraction
    
    def rotation_curve_visible(self, r: float) -> float:
        """
        可见物质的旋转曲线
        
        V(r) = √(GM(r) / r)
        
        对于指数盘: M(r) ~ (1 - e^(-r/R_d)) r²
        
        Args:
            r: 距中心距离 (kpc)
        
        Returns:
            旋转速度 (km/s)
        """
        G = AstroConstants.G
        r_m = r * 1e3 * AstroConstants.PARSEC
        
        # 指数盘模型
        R_d = self.R / 3  # 标度长度
        M_r = self.M_visible * (1 - np.exp(-r_m / R_d)) * (r_m / self.R)**2
        
        V = np.sqrt(G * M_r / r_m)
        
        return V / 1e3  # m/s → km/s
    
    def rotation_curve_total(self, r: float) -> float:
        """
        总旋转曲线（包括暗物质）
        
        V²(r) = V²_visible(r) + V²_dark(r)
        
        暗物质晕：NFW 轮廓
        
        Args:
            r: 距中心距离 (kpc)
        
        Returns:
            旋转速度 (km/s)
        """
        V_vis = self.rotation_curve_visible(r)
        
        # 暗物质贡献（简化 NFW）
        r_s = self.R / 10  # 特征半径 (kpc)
        
        if r < 0.1:
            V_dm = 0
        else:
            V_dm_sq = G * self.M_dark / (r * 1e3 * AstroConstants.PARSEC) * \
                     np.log(1 + r / r_s)
            V_dm = np.sqrt(V_dm_sq) / 1e3
        
        V_total = np.sqrt(V_vis**2 + V_dm**2)
        
        return V_total
    
    def virial_mass(self, velocity_dispersion: float) -> float:
        """
        维里质量
        
        M_vir = σ² R / G
        
        基于 ±0: 动能 ↔ 引力势能
        
        Args:
            velocity_dispersion: 速度弥散 (km/s)
        
        Returns:
            维里质量 (M_☉)
        """
        G = AstroConstants.G
        M_sun = AstroConstants.M_SUN
        
        sigma = velocity_dispersion * 1e3  # km/s → m/s
        
        M_vir = sigma**2 * self.R / G
        
        return M_vir / M_sun
    
    def tully_fisher_relation(self) -> float:
        """
        Tully-Fisher 关系
        
        L ~ V^α
        
        α ≈ 4 (旋涡星系)
        
        光度-旋转速度关系
        """
        V_max = self.rotation_curve_total(10)  # 10 kpc 处
        
        alpha = 4.0
        L_sun = AstroConstants.L_SUN
        
        # 归一化（银河系）
        V_MW = 220  # km/s
        L_MW = 2e10 * L_sun
        
        L = L_MW * (V_max / V_MW)**alpha
        
        return L / L_sun  # L_☉
    
    @staticmethod
    def milky_way() -> Galaxy:
        """
        银河系参数
        """
        return Galaxy(
            galaxy_type=Galaxy.GalaxyType.BARRED_SPIRAL,
            mass=1.5e12,      # M_☉
            radius=15,        # kpc
            dark_matter_fraction=0.90
        )


class DarkMatter:
    """
    暗物质
    
    不发光但有引力的物质
    
    基于 ±0 的理解:
    - 可见 ↔ 不可见
    - 重子物质 ↔ 非重子物质
    
    证据:
    - 星系旋转曲线平坦
    - 引力透镜
    - 宇宙微波背景
    """
    
    @staticmethod
    def rotation_curve_problem() -> str:
        """
        旋转曲线问题
        
        观测: V(r) ~ constant (平坦)
        预期: V(r) ~ 1/√r (Kepler)
        
        需要暗物质晕
        """
        return """
        旋转曲线问题:
        
        仅考虑可见物质:
            V(r) ∝ 1/√r  (Kepler 下降)
        
        实际观测:
            V(r) ≈ const  (平坦)
        
        解释:
            需要暗物质晕 ρ_DM(r) ∝ 1/r²
        """
    
    @staticmethod
    def nfw_profile(r: float, r_s: float, rho_0: float) -> float:
        """
        NFW 轮廓（暗物质密度分布）
        
        ρ(r) = ρ₀ / [(r/r_s)(1 + r/r_s)²]
        
        Navarro-Frenk-White profile
        
        Args:
            r: 距中心距离 (kpc)
            r_s: 特征半径 (kpc)
            rho_0: 归一化密度 (M_☉/kpc³)
        
        Returns:
            暗物质密度 (M_☉/kpc³)
        """
        x = r / r_s
        
        rho = rho_0 / (x * (1 + x)**2)
        
        return rho
    
    @staticmethod
    def candidates() -> List[str]:
        """
        暗物质候选
        """
        return [
            "WIMPs (Weakly Interacting Massive Particles)",
            "Axions",
            "Sterile neutrinos",
            "Primordial black holes",
            "MACHOs (Massive Compact Halo Objects)"
        ]


# ============================================================================
# Level 3: 宇宙学
# ============================================================================

class Cosmology:
    """
    宇宙学
    
    研究宇宙的起源、演化、结构
    
    基于 ±0 的理解:
    - 膨胀 ↔ 引力（宇宙动力学）
    - 物质 ↔ 能量（密度演化）
    - 均匀 ↔ 结构（扰动增长）
    """
    
    def __init__(self,
                 H0: float = AstroConstants.H0,
                 Omega_m: float = AstroConstants.OMEGA_M,
                 Omega_Lambda: float = AstroConstants.OMEGA_LAMBDA,
                 Omega_r: float = AstroConstants.OMEGA_R):
        """
        ΛCDM 宇宙学模型
        
        Args:
            H0: Hubble 常数 (km/s/Mpc)
            Omega_m: 物质密度参数
            Omega_Lambda: 暗能量密度参数
            Omega_r: 辐射密度参数
        """
        self.H0 = H0
        self.Omega_m = Omega_m
        self.Omega_Lambda = Omega_Lambda
        self.Omega_r = Omega_r
        self.Omega_k = 1 - Omega_m - Omega_Lambda - Omega_r
    
    def hubble_parameter(self, z: float) -> float:
        """
        Hubble 参数演化
        
        H(z) = H₀ √[Ω_m(1+z)³ + Ω_Λ + Ω_r(1+z)⁴ + Ω_k(1+z)²]
        
        Args:
            z: 红移
        
        Returns:
            H(z) (km/s/Mpc)
        """
        term = (self.Omega_m * (1 + z)**3 +
                self.Omega_Lambda +
                self.Omega_r * (1 + z)**4 +
                self.Omega_k * (1 + z)**2)
        
        H_z = self.H0 * np.sqrt(term)
        
        return H_z
    
    def age_of_universe(self) -> float:
        """
        宇宙年龄
        
        t₀ = ∫₀^∞ dz / [(1+z) H(z)]
        
        Returns:
            年龄 (年)
        """
        def integrand(z):
            return 1 / ((1 + z) * self.hubble_parameter(z))
        
        # 转换单位
        H0_SI = self.H0 * 1e3 / (1e6 * AstroConstants.PARSEC)  # s^-1
        
        t_H = 1 / H0_SI  # Hubble 时间 (s)
        
        integral, _ = quad(integrand, 0, np.inf)
        
        t_0 = t_H * integral
        
        # 秒 → 年
        t_0_years = t_0 / (365.25 * 24 * 3600)
        
        return t_0_years
    
    def critical_density(self, z: float = 0) -> float:
        """
        临界密度
        
        ρ_c(z) = 3H²(z) / (8πG)
        
        基于 ±0: 膨胀动能 ↔ 引力势能
        
        Args:
            z: 红移
        
        Returns:
            临界密度 (kg/m³)
        """
        G = AstroConstants.G
        H_z = self.hubble_parameter(z)
        
        # 转换单位
        H_SI = H_z * 1e3 / (1e6 * AstroConstants.PARSEC)  # s^-1
        
        rho_c = 3 * H_SI**2 / (8 * np.pi * G)
        
        return rho_c
    
    def comoving_distance(self, z: float) -> float:
        """
        共动距离
        
        D_C = c/H₀ ∫₀^z dz' / E(z')
        
        其中 E(z) = H(z)/H₀
        
        Args:
            z: 红移
        
        Returns:
            共动距离 (Mpc)
        """
        c = AstroConstants.C
        
        def integrand(z_prime):
            return 1 / self.hubble_parameter(z_prime) * self.H0
        
        integral, _ = quad(integrand, 0, z)
        
        D_C = (c / 1e3) * integral  # km/s * s/km = Mpc
        
        return D_C
    
    def luminosity_distance(self, z: float) -> float:
        """
        光度距离
        
        D_L = (1+z) D_C
        
        用于标准烛光（Ia 型超新星）
        
        Args:
            z: 红移
        
        Returns:
            光度距离 (Mpc)
        """
        D_C = self.comoving_distance(z)
        
        return (1 + z) * D_C
    
    def friedmann_equation(self, a: float) -> float:
        """
        Friedmann 方程
        
        (ȧ/a)² = (8πG/3) ρ - k/a² + Λ/3
        
        描述宇宙膨胀动力学
        
        基于 ±0: 膨胀 ↔ 引力减速 ↔ 暗能量加速
        
        Args:
            a: 标度因子
        
        Returns:
            (ȧ/a)² (H² 的归一化)
        """
        z = 1/a - 1
        H_a = self.hubble_parameter(z)
        
        return (H_a / self.H0)**2
    
    def acceleration_equation(self, a: float) -> float:
        """
        加速度方程
        
        ä/a = -(4πG/3)(ρ + 3P) + Λ/3
        
        基于 ±0: 物质减速(-) ↔ 暗能量加速(+)
        
        Args:
            a: 标度因子
        
        Returns:
            ä/a (归一化)
        """
        z = 1/a - 1
        
        # 简化：P_m ≈ 0, P_Λ = -ρ_Λ
        term = (-self.Omega_m * (1 + z)**3 / 2 +
                self.Omega_Lambda)
        
        return self.H0**2 * term
    
    def is_accelerating(self, z: float = 0) -> bool:
        """
        判断宇宙是否加速膨胀
        
        ä > 0 当 Ω_Λ > Ω_m/2
        """
        a = 1 / (1 + z)
        return self.acceleration_equation(a) > 0


class BigBang:
    """
    大爆炸理论
    
    宇宙从极热极密的奇点开始膨胀
    
    基于 ±0 的理解:
    - 奇点 ↔ 无限膨胀
    - 高温 ↔ 低温
    - 均匀 ↔ 结构
    """
    
    @staticmethod
    def timeline() -> Dict[float, str]:
        """
        宇宙演化时间线
        
        Returns:
            {时间(s): 事件}
        """
        return {
            0: "大爆炸奇点 (T → ∞)",
            1e-43: "Planck 时期 (量子引力)",
            1e-36: "暴胀开始",
            1e-32: "暴胀结束，再加热",
            1e-12: "电弱对称破缺",
            1e-6: "夸克-胶子等离子体 → 强子",
            1: "中微子退耦",
            180: "核合成开始 (BBN)",
            1200: "核合成结束",
            3.8e5 * 365.25 * 24 * 3600: "复合时期，CMB 形成 (z~1100)",
            1e8 * 365.25 * 24 * 3600: "第一代恒星形成",
            1e9 * 365.25 * 24 * 3600: "第一代星系形成",
            13.8e9 * 365.25 * 24 * 3600: "今天"
        }
    
    @staticmethod
    def nucleosynthesis() -> Dict[str, float]:
        """
        大爆炸核合成（BBN）
        
        t ~ 3-20 分钟
        T ~ 10⁹ K
        
        产生轻元素丰度
        
        Returns:
            质量比例
        """
        return {
            'H': 0.75,      # 氢
            'He-4': 0.25,   # 氦-4
            'D': 2e-5,      # 氘
            'He-3': 1e-5,   # 氦-3
            'Li-7': 1e-9    # 锂-7
        }
    
    @staticmethod
    def cmb_temperature(z: float) -> float:
        """
        宇宙微波背景温度演化
        
        T(z) = T₀ (1 + z)
        
        Args:
            z: 红移
        
        Returns:
            温度 (K)
        """
        T_0 = AstroConstants.T_CMB
        
        return T_0 * (1 + z)
    
    @staticmethod
    def recombination_redshift() -> float:
        """
        复合红移
        
        z_rec ~ 1100
        
        电子与质子结合成氢原子
        宇宙从不透明 → 透明
        """
        return 1100
    
    @staticmethod
    def matter_radiation_equality() -> float:
        """
        物质-辐射相等红移
        
        z_eq ~ 3400
        
        ρ_m = ρ_r
        
        从辐射主导 → 物质主导
        """
        Omega_m = AstroConstants.OMEGA_M
        Omega_r = AstroConstants.OMEGA_R
        
        z_eq = Omega_m / Omega_r - 1
        
        return z_eq


class Inflation:
    """
    暴胀理论
    
    宇宙极早期的指数膨胀
    
    解决:
    - 平坦性问题
    - 视界问题
    - 磁单极问题
    
    基于 ±0 的理解:
    - 标量场真空能量驱动膨胀
    - 慢滚 ↔ 快速衰减
    """
    
    def __init__(self, phi_0: float = 1.0, V_0: float = 1.0):
        """
        单场暴胀模型
        
        Args:
            phi_0: 标量场初值 (Planck 单位)
            V_0: 势能标度 (Planck 单位)
        """
        self.phi_0 = phi_0
        self.V_0 = V_0
    
    def potential(self, phi: float) -> float:
        """
        暴胀势能
        
        V(φ) = V₀ φ²  (简单二次势)
        
        Args:
            phi: 标量场值
        
        Returns:
            势能
        """
        return self.V_0 * phi**2 / 2
    
    def slow_roll_parameters(self, phi: float) -> Tuple[float, float]:
        """
        慢滚参数
        
        ε = (1/2) (V'/V)²
        η = V''/V
        
        慢滚条件: ε, |η| << 1
        
        Args:
            phi: 标量场值
        
        Returns:
            (ε, η)
        """
        V = self.potential(phi)
        V_prime = self.V_0 * phi
        V_double_prime = self.V_0
        
        epsilon = 0.5 * (V_prime / V)**2
        eta = V_double_prime / V
        
        return (epsilon, eta)
    
    def number_of_e_folds(self, phi_start: float, phi_end: float) -> float:
        """
        e-折叠数
        
        N = ∫ H dt = ∫ (V/V') dφ
        
        测量暴胀持续时间
        
        Args:
            phi_start: 起始场值
            phi_end: 结束场值
        
        Returns:
            N
        """
        # 对于二次势
        N = (phi_start**2 - phi_end**2) / (2 * self.V_0)
        
        return N
    
    @staticmethod
    def solve_horizon_problem() -> str:
        """
        解决视界问题
        """
        return """
        视界问题:
        
        问题: CMB 各向同性，但因果联系区域太小
        
        解决: 暴胀使当前视界内的区域来自
              暴胀前的一个小因果区域
              
              指数膨胀: a(t) ~ e^(Ht)
              N ~ 60 e-folds 足够
        """
    
    @staticmethod
    def primordial_fluctuations() -> str:
        """
        原初扰动
        
        量子涨落 → 经典扰动 → 结构种子
        """
        return """
        原初扰动:
        
        1. 暴胀期间量子涨落
        2. 膨胀冻结为经典扰动
        3. 功率谱: P(k) ~ k^(n_s-1)
           标量谱指数: n_s ≈ 0.96
        4. 成为结构形成的种子
        """


class StructureFormation:
    """
    宇宙结构形成
    
    从微小密度扰动到星系、星系团
    
    基于 ±0 的理解:
    - 引力不稳定性（扰动增长）
    - 过密区(+) ↔ 欠密区(-)
    """
    
    @staticmethod
    def jeans_length(rho: float, T: float, mu: float = 0.6) -> float:
        """
        Jeans 长度
        
        λ_J = √(π c_s² / G ρ)
        
        大于此尺度的扰动会坍缩
        
        基于 ±0: 压力支撑(+) ↔ 引力(-)
        
        Args:
            rho: 密度 (kg/m³)
            T: 温度 (K)
            mu: 平均分子量
        
        Returns:
            Jeans 长度 (m)
        """
        G = AstroConstants.G
        k_B = AstroConstants.K_B
        m_H = 1.67e-27  # kg
        
        # 声速
        c_s = np.sqrt(k_B * T / (mu * m_H))
        
        lambda_J = np.sqrt(np.pi * c_s**2 / (G * rho))
        
        return lambda_J
    
    @staticmethod
    def growth_factor(z: float, Omega_m: float = 0.315) -> float:
        """
        线性增长因子
        
        D(z) / D(0) ∝ (1+z)^(-1) (物质主导)
        
        扰动振幅增长
        
        Args:
            z: 红移
            Omega_m: 物质密度参数
        
        Returns:
            D(z)/D(0)
        """
        # 简化（Einstein-de Sitter）
        D_ratio = 1 / (1 + z)
        
        return D_ratio
    
    @staticmethod
    def power_spectrum(k: float, A_s: float = 2.1e-9, n_s: float = 0.96) -> float:
        """
        功率谱
        
        P(k) = A_s k^(n_s - 1)
        
        描述密度扰动的统计性质
        
        Args:
            k: 波数 (h/Mpc)
            A_s: 振幅
            n_s: 标量谱指数
        
        Returns:
            P(k)
        """
        return A_s * k**(n_s - 1)
    
    @staticmethod
    def press_schechter_mass_function(M: float, z: float = 0) -> float:
        """
        Press-Schechter 质量函数
        
        dn/dM: 质量为 M 的晕的数密度
        
        Args:
            M: 晕质量 (M_☉)
            z: 红移
        
        Returns:
            dn/dM (Mpc^-3 M_☉^-1)
        """
        # 简化公式
        delta_c = 1.686  # 临界过密度
        sigma_M = 1.0    # 归一化（应从功率谱计算）
        
        rho_m = AstroConstants.OMEGA_M * 2.78e11  # M_☉ / h Mpc³
        
        dn_dM = (rho_m / M**2) * np.exp(-delta_c**2 / (2 * sigma_M**2))
        
        return dn_dM


# ============================================================================
# Level 4: 黑洞
# ============================================================================

class BlackHole:
    """
    黑洞
    
    时空极度弯曲，光无法逃逸
    
    基于 ±0 的理解:
    - 视界内 ↔ 视界外（信息边界）
    - Hawking 辐射：黑洞蒸发
    """
    
    def __init__(self, mass: float):
        """
        Schwarzschild 黑洞
        
        Args:
            mass: 质量 (M_☉)
        """
        self.M = mass * AstroConstants.M_SUN
    
    def schwarzschild_radius(self) -> float:
        """
        Schwarzschild 半径（事件视界）
        
        R_s = 2GM/c²
        """
        G = AstroConstants.G
        c = AstroConstants.C
        
        return 2 * G * self.M / c**2
    
    def photon_sphere_radius(self) -> float:
        """
        光子球半径
        
        R_ph = 3GM/c²
        
        光的不稳定圆轨道
        """
        return 1.5 * self.schwarzschild_radius()
    
    def innermost_stable_circular_orbit(self) -> float:
        """
        最内稳定圆轨道（ISCO）
        
        R_ISCO = 6GM/c²
        """
        return 3 * self.schwarzschild_radius()
    
    def hawking_temperature(self) -> float:
        """
        Hawking 温度
        
        T_H = ℏc³ / (8πGMk_B)
        
        黑洞辐射温度
        
        基于 ±0: 视界附近的虚粒子对分离
        """
        G = AstroConstants.G
        c = AstroConstants.C
        hbar = AstroConstants.H_BAR
        k_B = AstroConstants.K_B
        
        T_H = hbar * c**3 / (8 * np.pi * G * self.M * k_B)
        
        return T_H
    
    def hawking_luminosity(self) -> float:
        """
        Hawking 光度
        
        L = (ℏc⁶) / (15360πG²M²)
        """
        G = AstroConstants.G
        c = AstroConstants.C
        hbar = AstroConstants.H_BAR
        
        L = hbar * c**6 / (15360 * np.pi * G**2 * self.M**2)
        
        return L
    
    def evaporation_time(self) -> float:
        """
        蒸发时间
        
        τ ~ M³ / (ℏc⁴)
        
        黑洞完全蒸发所需时间
        """
        G = AstroConstants.G
        c = AstroConstants.C
        hbar = AstroConstants.H_BAR
        
        tau = 5120 * np.pi * G**2 * self.M**3 / (hbar * c**4)
        
        return tau / (365.25 * 24 * 3600)  # 秒 → 年
    
    def entropy(self) -> float:
        """
        Bekenstein-Hawking 熵
        
        S_BH = (k_B c³ / 4ℏG) A
        
        其中 A = 4π R_s² 是视界面积
        
        基于 ±0: 熵 ∝ 面积（全息原理）
        """
        k_B = AstroConstants.K_B
        c = AstroConstants.C
        G = AstroConstants.G
        hbar = AstroConstants.H_BAR
        
        R_s = self.schwarzschild_radius()
        A = 4 * np.pi * R_s**2
        
        S_BH = (k_B * c**3 / (4 * hbar * G)) * A
        
        return S_BH / k_B  # 单位 k_B


class GravitationalWave:
    """
    引力波
    
    时空的涟漪
    
    基于 ±0 的理解:
    - 拉伸 ↔ 压缩（应变）
    - 质量四极矩变化
    """
    
    @staticmethod
    def chirp_mass(m1: float, m2: float) -> float:
        """
        啁啾质量
        
        ℳ = (m₁m₂)^(3/5) / (m₁+m₂)^(1/5)
        
        决定引力波频率演化
        
        Args:
            m1, m2: 双星质量 (M_☉)
        
        Returns:
            啁啾质量 (M_☉)
        """
        M_chirp = (m1 * m2)**(3/5) / (m1 + m2)**(1/5)
        
        return M_chirp
    
    @staticmethod
    def orbital_frequency(m1: float, m2: float, separation: float) -> float:
        """
        轨道频率
        
        f_orb = (1/2π) √[G(m₁+m₂)/a³]
        
        Args:
            m1, m2: 质量 (M_☉)
            separation: 轨道分离 (R_☉)
        
        Returns:
            轨道频率 (Hz)
        """
        G = AstroConstants.G
        M_sun = AstroConstants.M_SUN
        R_sun = AstroConstants.R_SUN
        
        M_total = (m1 + m2) * M_sun
        a = separation * R_sun
        
        f_orb = (1 / (2 * np.pi)) * np.sqrt(G * M_total / a**3)
        
        return f_orb
    
    @staticmethod
    def strain_amplitude(M_chirp: float, distance: float, frequency: float) -> float:
        """
        应变振幅
        
        h ~ (Gℳc/c²) / d × (πGℳf/c³)^(2/3)
        
        Args:
            M_chirp: 啁啾质量 (M_☉)
            distance: 距离 (Mpc)
            frequency: 频率 (Hz)
        
        Returns:
            应变 h (无量纲)
        """
        G = AstroConstants.G
        c = AstroConstants.C
        M_sun = AstroConstants.M_SUN
        
        M_c = M_chirp * M_sun
        d = distance * 1e6 * AstroConstants.PARSEC
        
        h = (G * M_c / c**2) / d * (np.pi * G * M_c * frequency / c**3)**(2/3)
        
        return h
    
    @staticmethod
    def merger_time(m1: float, m2: float, separation: float) -> float:
        """
        并合时间
        
        τ = (5c⁵/256G³) × a⁴/(m₁m₂(m₁+m₂))
        
        Args:
            m1, m2: 质量 (M_☉)
            separation: 初始分离 (R_☉)
        
        Returns:
            并合时间 (年)
        """
        G = AstroConstants.G
        c = AstroConstants.C
        M_sun = AstroConstants.M_SUN
        R_sun = AstroConstants.R_SUN
        
        m1_kg = m1 * M_sun
        m2_kg = m2 * M_sun
        a = separation * R_sun
        
        tau = (5 * c**5 / (256 * G**3)) * \
              (a**4 / (m1_kg * m2_kg * (m1_kg + m2_kg)))
        
        return tau / (365.25 * 24 * 3600)  # 秒 → 年


# ============================================================================
# 综合演示系统
# ============================================================================

class AstrophysicsDemo:
    """天体物理综合演示"""
    
    @staticmethod
    def demonstrate_stellar_structure():
        """演示恒星结构"""
        print("\n" + "=" * 70)
        print("恒星结构")
        print("=" * 70)
        
        # 太阳
        sun = StellarStructure(
            mass=1.0,
            radius=1.0,
            luminosity=1.0,
            composition={'H': 0.70, 'He': 0.28, 'Z': 0.02}
        )
        
        print(f"\n太阳:")
        print(f"  质量: {sun.M / AstroConstants.M_SUN:.2f} M_☉")
        print(f"  半径: {sun.R / AstroConstants.R_SUN:.2f} R_☉")
        print(f"  光度: {sun.L / AstroConstants.L_SUN:.2f} L_☉")
        
        T_c = sun.central_temperature()
        T_surf = sun.surface_temperature()
        
        print(f"\n温度:")
        print(f"  中心温度: {T_c:.2e} K")
        print(f"  表面温度: {T_surf:.0f} K")
        
        # 求解结构
        structure = sun.solve_structure(num_points=50)
        
        print(f"\n结构剖面:")
        print(f"  中心密度: {structure['rho'][0]:.2e} kg/m³")
        print(f"  中心压力: {structure['P'][0]:.2e} Pa")
        
        # 核反应
        pp = NuclearFusion.pp_chain()
        cno = NuclearFusion.cno_cycle()
        
        print(f"\n核反应:")
        print(f"  PP 链: {pp['reaction']}")
        print(f"    能量释放: {pp['Q_value']:.1f} MeV")
        print(f"    温度范围: {pp['temperature_range']}")
        
        print(f"\n  CNO 循环: {cno['reaction']}")
        print(f"    能量释放: {cno['Q_value']:.1f} MeV")
        print(f"    主导质量: {cno['dominant_in']}")
    
    @staticmethod
    def demonstrate_stellar_evolution():
        """演示恒星演化"""
        print("\n" + "=" * 70)
        print("恒星演化")
        print("=" * 70)
        
        masses = [0.5, 1.0, 10.0, 30.0]
        
        print(f"\n主序寿命:")
        for M in masses:
            tau = StellarEvolution.main_sequence_lifetime(M)
            L = StellarEvolution.mass_luminosity_relation(M)
            final = StellarEvolution.final_state(M)
            
            print(f"  M = {M:.1f} M_☉:")
            print(f"    光度: {L:.2e} L_☉")
            print(f"    寿命: {tau:.2e} 年")
            print(f"    最终归宿: {final.value}")
        
        # 特殊质量限制
        M_Ch = StellarEvolution.chandrasekhar_limit()
        M_TOV = StellarEvolution.tolman_oppenheimer_volkoff_limit()
        
        print(f"\n质量极限:")
        print(f"  Chandrasekhar 极限: {M_Ch:.2f} M_☉")
        print(f"    (白矮星最大质量)")
        print(f"  TOV 极限: {M_TOV:.2f} M_☉")
        print(f"    (中子星最大质量)")
        
        # Schwarzschild 半径
        print(f"\nSchwarzschild 半径:")
        for M in [1, 10, 100]:
            R_s = StellarEvolution.schwarzschild_radius(M)
            print(f"  M = {M} M_☉: R_s = {R_s/1e3:.2f} km")
    
    @staticmethod
    def demonstrate_galaxies():
        """演示星系"""
        print("\n" + "=" * 70)
        print("星系与暗物质")
        print("=" * 70)
        
        # 银河系
        milky_way = Galaxy.milky_way()
        
        print(f"\n银河系:")
        print(f"  类型: {milky_way.type.value}")
        print(f"  总质量: {milky_way.M_total / AstroConstants.M_SUN:.2e} M_☉")
        print(f"  可见物质: {milky_way.M_visible / AstroConstants.M_SUN:.2e} M_☉")
        print(f"  暗物质: {milky_way.M_dark / AstroConstants.M_SUN:.2e} M_☉")
        print(f"  暗物质比例: {milky_way.f_dm:.2%}")
        
        # 旋转曲线
        print(f"\n旋转曲线:")
        r_vals = [1, 5, 10, 20]  # kpc
        for r in r_vals:
            V_vis = milky_way.rotation_curve_visible(r)
            V_total = milky_way.rotation_curve_total(r)
            
            print(f"  r = {r:2d} kpc:")
            print(f"    可见物质: V = {V_vis:.0f} km/s")
            print(f"    总（含暗物质）: V = {V_total:.0f} km/s")
        
        # 暗物质证据
        print(f"\n暗物质证据:")
        print(DarkMatter.rotation_curve_problem())
        
        print(f"\n暗物质候选:")
        for candidate in DarkMatter.candidates():
            print(f"  • {candidate}")
    
    @staticmethod
    def demonstrate_cosmology():
        """演示宇宙学"""
        print("\n" + "=" * 70)
        print("宇宙学")
        print("=" * 70)
        
        # ΛCDM 模型
        cosmo = Cosmology()
        
        print(f"\nΛCDM 宇宙学参数:")
        print(f"  Hubble 常数: H₀ = {cosmo.H0:.1f} km/s/Mpc")
        print(f"  物质密度: Ω_m = {cosmo.Omega_m:.3f}")
        print(f"  暗能量密度: Ω_Λ = {cosmo.Omega_Lambda:.3f}")
        print(f"  辐射密度: Ω_r = {cosmo.Omega_r:.2e}")
        print(f"  曲率: Ω_k = {cosmo.Omega_k:.3f}")
        
        # 宇宙年龄
        t_0 = cosmo.age_of_universe()
        print(f"\n宇宙年龄: t₀ = {t_0/1e9:.2f} Gyr")
        
        # 临界密度
        rho_c = cosmo.critical_density(0)
        print(f"临界密度: ρ_c = {rho_c:.2e} kg/m³")
        
        # Hubble 参数演化
        print(f"\nHubble 参数演化:")
        z_vals = [0, 0.5, 1, 2, 5, 10]
        for z in z_vals:
            H_z = cosmo.hubble_parameter(z)
            print(f"  z = {z:2d}: H(z) = {H_z:.1f} km/s/Mpc")
        
        # 距离
        print(f"\n距离测度:")
        z_test = [0.1, 0.5, 1.0, 2.0]
        for z in z_test:
            D_C = cosmo.comoving_distance(z)
            D_L = cosmo.luminosity_distance(z)
            print(f"  z = {z:.1f}:")
            print(f"    共动距离: D_C = {D_C:.0f} Mpc")
            print(f"    光度距离: D_L = {D_L:.0f} Mpc")
        
        # 加速膨胀
        is_accel = cosmo.is_accelerating(0)
        print(f"\n当前宇宙加速膨胀: {is_accel}")
    
    @staticmethod
    def demonstrate_big_bang():
        """演示大爆炸"""
        print("\n" + "=" * 70)
        print("大爆炸理论")
        print("=" * 70)
        
        # 时间线
        timeline = BigBang.timeline()
        print(f"\n宇宙演化时间线:")
        for t, event in list(timeline.items())[:8]:
            if t < 1:
                print(f"  t = {t:.2e} s: {event}")
            elif t < 3600:
                print(f"  t = {t:.0f} s: {event}")
            else:
                print(f"  t = {t/(365.25*24*3600):.2e} yr: {event}")
        
        # 核合成
        bbn = BigBang.nucleosynthesis()
        print(f"\n大爆炸核合成（t ~ 3 min）:")
        for element, abundance in bbn.items():
            print(f"  {element:6s}: {abundance:.2e}")
        
        # CMB
        z_rec = BigBang.recombination_redshift()
        T_rec = BigBang.cmb_temperature(z_rec)
        
        print(f"\n宇宙微波背景:")
        print(f"  复合红移: z_rec = {z_rec:.0f}")
        print(f"  复合温度: T_rec = {T_rec:.0f} K")
        print(f"  今天温度: T_CMB = {AstroConstants.T_CMB:.3f} K")
        
        # 物质-辐射相等
        z_eq = BigBang.matter_radiation_equality()
        print(f"\n物质-辐射相等:")
        print(f"  z_eq = {z_eq:.0f}")
    
    @staticmethod
    def demonstrate_black_holes():
        """演示黑洞"""
        print("\n" + "=" * 70)
        print("黑洞物理")
        print("=" * 70)
        
        masses = [1, 10, 1e6, 1e9]  # M_☉
        
        for M in masses:
            bh = BlackHole(M)
            
            R_s = bh.schwarzschild_radius()
            R_isco = bh.innermost_stable_circular_orbit()
            T_H = bh.hawking_temperature()
            tau_evap = bh.evaporation_time()
            S = bh.entropy()
            
            print(f"\nM = {M:.0e} M_☉:")
            print(f"  Schwarzschild 半径: R_s = {R_s/1e3:.2e} km")
            print(f"  ISCO: R_ISCO = {R_isco/1e3:.2e} km")
            print(f"  Hawking 温度: T_H = {T_H:.2e} K")
            print(f"  蒸发时间: τ = {tau_evap:.2e} 年")
            print(f"  熵: S = {S:.2e} k_B")
        
        # 引力波
        print(f"\n引力波:")
        m1, m2 = 30, 30  # M_☉
        M_chirp = GravitationalWave.chirp_mass(m1, m2)
        
        print(f"  双黑洞系统: {m1} M_☉ + {m2} M_☉")
        print(f"  啁啾质量: ℳ = {M_chirp:.2f} M_☉")
        
        separation = 1000  # R_☉
        f_orb = GravitationalWave.orbital_frequency(m1, m2, separation)
        tau_merge = GravitationalWave.merger_time(m1, m2, separation)
        
        print(f"  初始分离: a = {separation} R_☉")
        print(f"  轨道频率: f = {f_orb:.2e} Hz")
        print(f"  并合时间: τ = {tau_merge:.2e} 年")


def main():
    """主函数：完整演示"""
    
    print("=" * 70)
    print("DaoMath Astrophysics")
    print("从 ±0 构造宇宙理论")
    print("=" * 70)
    
    # 1. 恒星结构
    AstrophysicsDemo.demonstrate_stellar_structure()
    
    # 2. 恒星演化
    AstrophysicsDemo.demonstrate_stellar_evolution()
    
    # 3. 星系与暗物质
    AstrophysicsDemo.demonstrate_galaxies()
    
    # 4. 宇宙学
    AstrophysicsDemo.demonstrate_cosmology()
    
    # 5. 大爆炸
    AstrophysicsDemo.demonstrate_big_bang()
    
    # 6. 黑洞
    AstrophysicsDemo.demonstrate_black_holes()
    
    print("\n" + "=" * 70)
    print("天体物理构造完成！")
    print("=" * 70)


if __name__ == "__main__":
    main()
